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Abstract. - We have built a database of ab-initio total energies for elemental Boron in over 60 
hypothetical crystal structures of varying coordination Z, such that every atom is equivalent. 
Fitting to each subset with a particular Z, we extract a classical effective potential, written as 
a sum over coordination shells and dominated by three-atom (bond angle dependent) terms. In 
the case Z — 6 (lowest in energy and most relevant), the classical potential has a typical error 
of 0.1 eV/atom, and favours the "inverted-umbrells" environment seen in real Boron. 



Introduction. - Boron stands out among the elements for the complexity and polymorphous 
variety of its structures fl. The structures are networks built from icosahedrally symmetric 
clusters, and the typical "inverted-umbrella" coordination shell is asymmetrical, containing 
five neighbors on one side and one on the other. The stable phase (/3-B) is currently modeled 
with 320 atoms/unit cell. Of the metastable allomorphs jL|, only a±2 and T50 have solved 
structures the latter of which is stabilized by 4% N atoms jl. The simplest phase 
a\2 contains 3-centered bonds both within and between icosahedra || . It was thought that 
bonding within icosahedra is metallic and weaker, while inter-icosahedral bonds are covalent 
and strong @, |6); however, recent experiments have questioned this behavior |t], @|. 

Boron is the only plausible candidate for a single-element or covalent quasicrystal, which 
can be speculatively modeled H |hJ on either ai2-B or /3-B. This idea inspired the discovery 
of a new Boron phase |ll| , as well as the prediction of Boron nanotubes [Q . 

Our motivation is to compare the energies of alternate Boron structures. Selected crystal 



structures |i| |13J, finite icosahedral clusters ||, microtube segments, or sheet fragments |12 
were compared by direct ab-initio calculations. Such computations, however, are limited to 
O(10 2 ) atoms whereas real Boron has more atoms per unit cell, and quasicrystal models require 
~ 10 5 atoms. Furthermore Monte Carlo and/or molecular dynamics simulations are desirable, 
especially for the liquid [O. 
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A tractable approximation is required of the ab-initio total energy as a function of atomic 
positions, such as a classical (many-atom) potential. This potential must reasonably represent 
the energy not only for the ground-state structure with slight distortions, but also for relatively 
high-energy local environments, such as occur in defects or sometimes in complex ground states, 
We know no previous attempts at such a potential. (The closest precursor of our work is Lee's 
study of a few Boron structures [fist , comparing atomic-orbital total energies with a simpler 
geometrical function of the closed circuits in the Boron network.) 

We start from a general form of potential into which all Silicon potentials jlfj, [IT], [l8|, 19 



could be cast, a sum E to tai — Yli Ei over interactions local to each coordination shell: The 
site energy Ei of atom i with coordination number Zi is broken into terms depending on one, 
two, etc. of its neighbors j, k, 

Ei =Y,f^) + H G ^ r ^nk,efl) + ... (l) 

(i) 

where Tij is the distance from atom i to j, and 6*) fe the angle formed by two neighbours j, A; 
and center i. (We shall assume later that j, k are restricted to "nearest" neighbours of site i.) 

Those Si potentials mostly attempted to capture the role of Z and the radial dependences. 
The dependence on bond angles (with fixed coordination and bond radius) is not so assured, 
even for Si. Three-body terms may be implicit, as each bond's effective Z depends on the 
surrounding atoms [|l7| Q. Bazant et al |l9| critiqued the angle dependence of prior potentials, 
but they too addressed it indirectly, by assuming various angular forms and checking which of 
these gave the most transferable radial behaviours. 

In this letter, we develop the opposite approach: we isolate the angular dependence without 
a priori assumptions of analytic form, by fitting to a large database in which non-angular vari- 
ables are held fixed. We have no hopes that a potential should be transferable to a structures 
with very different local order; thus our database must be a large family of structures, among 



which bond angles and coordination numbers vary independently. I 1 ) This is not the case 
for previous databases of periodic structures. The reason is, in part, that at Boron's typical 
coordination Z = 6, a much greater variety of coordination shells is plausible than at ordinary 
covalent (Z — 3-4) or metallic (Z w 12) coordinations. 

Consider a database of structures in which all bond lengths are R, all coordination numbers 
arc Z, and non-pair interactions are limited to nearest neighbors. Then eq. (|l]) reduces to 

Ei = J2 fzinj) + G z{R, © + h z (i) (2) 

where Gz [R, 9) = Gz (R, R, 9) , and hz gathers all interactions beyond the three-body (or 
rather two-neighbor) terms. It is not obvious a priori what form hz should have; we found 
empirically (see below) that four-body (three-neighbor) terms are not needed, but a Z-neighbor 
term hz(i) = cz{R)t;{i) is needed. It is proportional to the "asymmetry" £(i) of coordination 
shell i (originally introduced to characterize dangling bonds p8| ): 

where is the vector from atom i to its neighbour j. (The denominator is the mean 
nearest-neighbour distance.) 

( 1 ) We can quantify how well our database samples the space of possible coordination shells, by 
use of a metric which defines a distance in this space [^| . 
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Fig. 1. - Some less familiar Z& structures. Bravais lattice vertices are indicated by *. For layered 
structures, atoms • are nearer to the viewer than o. (a) Kagome, (b) i-dodec, (c) P[i, (d) Pa, (e) two 
layers in T-lattice (the full cubic structure is an ABC stacking of Kagome layers) (f) two layers in 
SC-co (full structure is • • o stacking), (g) SC-ico, (h) tube-<r: the a- lattice at top is rolled into tubes 
(identifying the points connected by the dashed curve) which are seen end-on at bottom. 



Now, restrict the database further so that every structure is "uniform" i.e. all sites are 
crystallographically equivalent and have the same local environment; We investigate those 
uniform structures in which all nearest neighbour distances are made equal (or nearly so), the 
better to separate the i?-dependence from the angular effects. 

Then Ei must be equated with the total energy per atom, as found from an LDA calculation. 
Now assume that, within the database, only a discrete set of inter-neighbor angles is possible, 
9 a (a = 1,2, ...); we realize this approximately by simply dividing the range of 9 into several 
bins. Then - if the structures sufficiently outnumber the angular bins - one obtains Gz{R, 9 a ) 
for each angle by a simple linear fit. For each value of the scale R, the coefficients {Gz(R, 9 a )} 
and Cz(R) satisfy a set of linear equations: 

E m (R) - Zf z (R) = N m (9 a )G z (R, Oa)+C Z (R)£,m (4) 
a 

Here m runs over the structures of the same Z, and N m {9 a ) is the number times angle 9 a 
occurs in the coordination shell of structure m. The whole procedure can be repeated for each 
Z. 



Structural database. - For our database, we adopted or invented over 60 uniform struc- 
tures ]2C| , spanning coordination numbers Z from 3 to 12, and exhibiting a variety of angular 
patterns for each Z; We shall present details only for 6-coordinated (Z6) structures, which 
are relevant to the real phases: 80% of the sites in /3-B, and 50% in ai2(B), have Z — 6. The 
same procedure was carried out for Z — 4, 5 and 7. ||(| 

Our Z6 structures, of course, include the simple cubic (SC) and the triangular (tri) lattices 
(parentheses give codes to be used henceforth). We systematically constructed more Z6 
structures by stacking Z4 planar lattices, such as the 3.4.6.4 pattern of inter-penetrating 
dodecagons (i-dodec) or the 3.6.3.6 lattice (Kagome), see Figure |l](a,b). We also stacked Z5 
planar lattices such as 3 2 .4.3.4 (called a) or the 3 3 .4 2 (called /i), puckering alternate layers so 
that each vertex gained one more neighbour from the layer above or below: this produced the 
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Table I. - Local environments of Z6 structures. 



Structure 


N 2 


N 3 


60° 


90° 


108° 


N(9 a ) 
120° 


135° 


150° 


180° 




SC 


12 


8 





12 














3 





tri 





G 


G 








G 








3 





SC-co 


7 


6 


2 


7 





2 


4 








0.586 


T-lattice 





12 


G 








G 








3 





i-dodec 


10 


6 


1 


10 





1 





2 


1 


0.732 


Kagome 


8 


4 


2 


8 





2 








3 





SC-ico 


1 


10 


5 


1 


7 








2 





1.854 


Pa 


4 b 


12 c 


3 


3 d 


4 d 


l d 


2 d 


2 d 





0.263 




5 b 


10 c 


3 


4 


3 d 


2 


2 d 





1 


0.771 


tube-tri 





G a 


G 








6 b 


1 





2 d 


1.163 


tube-cr a 


5 b 


7 C 


3 


5 b 


2 b 


2 





3 b 





0.936 



Z6 puckered (Per) and {P/i) structures, respectively as in Figure ||(c,d). As far as possible, 
the inter- and intra-unit bond lengths in our packings and stackings were made equal. 

Figure |l](e,f,g) show inherently three-dimensional networks: a simple cubic array of cuboc- 
tahedra (SC-co); the (T-lattice) - also called "pyrochlore lattice" - consisting of a diamond 



lattice's bond midpoints; and icosahedra placed on a simple cubic lattice (SC-ico). Q We also 
rolled a triangular lattice into a single infinite nanotubc |l^] , with a circumference of 8 edges 
(tube-tri, not illustrated). Finally, we rolled the Zh 3 2 .4.3.4 (a) lattice into a tube whose 
circumference is the dotted arc in Figure |l|(h), and packing such tubes in a square array. This 
made the "tube-cr" structure, in which one neighbor of each atom belongs to an adjacent tube. 
Table [| summarizes geometric data on the coordination environments of these Z6 structures. 



( 3 ) Here N2 and N3 are the number of neighbours at at distances ~ y/2R and ~ \/3R, 
respectively. The other columns directly determine our potential (Q): the number N(0 a ) in 
each coordination shell of "two- neighbour" angles 6 a , and the asymmetry £ of the coordination 
shell defined by (||). 

LDA Results. - We performed an ab initio total energy calculation for each structure in 
our database, in the local density approximation with extended norm and hardness conserving 
pseudo-potentials [^l] . We used all planewaves with kinetic energy up to 54.5 Ry. The Brillouin 
zone is sampled with a k-point density of at least (16.3A" 1 ) 3 . Band structure energies are 
converged to within 10~ 6 Ry in each calculation. The convergence with respect to k-point 
density shows a precision of 0.03 eV for energy comparison among different structures. 

We varied the lattice constant of each structure (a uniform scale factor, without any 
relaxation of positions). Figure ||(a) collects the cohesive energy (per atom) of each Z6 
structure as a function of the nearest neighbour distance R. The inset shows the single- well- like 
shape for R up to 5A for three representative structures, SC, tri, and SC-ico; beyond R = 3A, 
the curves for all structures converge, approaching the non-bonding limit. The main figure 

( 2 ) SC-ico was introduced in Ref. |l5[ ], but misidentified there as "primitive orthorhombic." 

( 3 ) In the table, footnote a means nearest (~ R) neighbours are not at identical distances; b means 
N2 has a tolerance of 0.14i? in distance; c means neighbors at various distances between y/2R and 
V3R are binned in N3; d means 9 a has a tolerance of 6° in angle, and the "180°" bin in tube-tri is 
actually at 169°. 



W.-J. ZHU et al: POTENTIALS FOR 6-COORDINATED BORON ETC. 



5 




Fig. 2. - (a). Cohesive energy Ei( R) for Z6 structures. In the inset, the three representative structures 
are SC (line), tri (dotted), and SC-ico (long-dash). The region showing the energy wells is magnified in 
the main figure. The structures are, in order of lowest to highest energy at R = 1.8A, ideal-ai2 (line), 
SC-ico (long-dash), SC-co (□), tube-cr (o), T-lattice (line), P-fi (long-dash), tri (short-dash), P-cr (x), 
tube-tri (line), SC (line), i-dodec (o), and Kagome (x). (b). Angular and asymmetry potentials. The 
fitted G(R,9) is plotted with symbols {*, o, x,+,o}, corresponding to R — {1.6, 1.7, 1.8, 1.9, 2.0} A; 
the first and last sets are connected by solid and dashed lines, respectively. The coefficient c(R) (see 
eqs. ([]) and is shown by the same symbols, on the same scale, in the "Asym" column. 



enlarges the physically relevant range R € [1.6, 2.0]A, to highlight the differences among the 
structures. The lowest energy curve shown is the "ideal-ai2" phase, topologically same as ai2 
but with all nearest neighbour distances set equal. (This was omitted from our Z6 database 
because it has two kinds of sites, of which one has Z — 7.) The Z6 uniform structure of 
lowest energy is SC-ico, another way of connecting B12 icosahedra so each atom has one 
inter-icosahedral bond. Slightly higher is SC-co, which can be viewed analogously as a packing 
of B12 cuboctahedra in place of icosahedra, with each atom having two inter-cluster bonds. 

Fitting of an effective potential. ~ The simplest effective potential would be a single-well 
two-body potential of typical shape, so the nearest-neighbour radius lies close to the minimum 
of the well, while further neighbours contribute at the tail of the well. Then structures of 
the same coordination would have similar energies (from the near-neighbor terms) with small 
differences due to the farther neighbors. This will not work in general [|l8|] and failed badly 
for our energy curves. Basically, one is asking the tail of fzif) in (@) to account for both the 
total energy when R — ► ^/2R or ^/ZR - a sizeable fraction of the well depth - as well as the 
dependence on N2 and N3 (see Table Q) among structures of the same Z - which is smaller, as 
seen in Figure @. We conclude that the two-body potential must be essentially truncated after 
the nearest-neighbor distance. It follows that the first term in (ph reduces to Zf(R); the energy 
differences among structures of the same coordination Z must be attributed to higher-order 
potentials. For each Z, we must choose a reference structure for which the higher-order terms 
are set to zero (thereby defining Zfz(R) to be its energy curve); for Z6 we chose the simple 
cubic (SC) structure. 

Thus we are fitting the E m (R) from Figure |2la) to the linear equations (||), in which 
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Tabic || (minus columns N2 and N3) constitutes the 11x8 matrix of coefficients. Eqs. (|J) are 
overdetermined (by three degrees of freedom); hence for each R, we can solve for Gs(R,9 a ) 
and ce(R) in a least-squares sense, without assuming any functional forms. The resulting fit 
(Figure||) snows the angular potential Gq(R,9) increases monotonically with 9. Both cq(R) 
and Gq(R, 9) show rapid decay as a function of R. Qualitatively similar (R, 9) dependences are 
found for Z = 4, 5, and 7 structures pel , except that for Z = 4 the sign of G4 is reversed. The 
monotonicity of the fitted Gz(R,9), for all four values of Z, argues for the physical validity 
of the result: a spurious fit should produce random fluctuations as a function of 9 a , since the 
value of Gz at each bin is an independent parameter in the fit. 



Discussion. - The first check on our potentials is that, among the 11 structures in our 
ZQ database, the total energies are in the right order (apart from exchange of certain close 
ones), with an average error of 0.14 eV/atom (for R in the bonding range). Furthermore the 
"inverted-umbrella" coordination shell of real Boron - which was not used in our database 
- was correctly found lower (by ~ leV) than any of the Z6 shells that were used. We also 
performed a Monte Carlo search of random Z6 environments, at the valid R = 1.7A and 
found the inverted-umbrella was the second-lowest in energy. F] Finally, our Gq disfavouring 
9 = 180° (see Fig. ||) correctly predicts that Boron in triangular sheets buckles, as found by 
ab-initio calculations [^2[ . 

Empirical tight-binding calculations are the successful intermediate between a fully ab-initio 
and our classical-potential approach [p2[ . It would be worthwhile to fit our database in this 
fashion, especially if longer-range, metallic interactions are found within the clusters, tubes, 
and sheets. Tight-binding (atomic orbital) total energies can also be expanded in moments 
related to circuits of three or four atoms. Structure selection and bond angles was studied in 
this way by Lee et al p5) for selected packings of B12 icosahedra. It would be interesting to 
analytically relate the above-mentioned circuits to the form of classical potentials found by us. 

Can our potentials, limited to fixed Z and R, be extended to arbitrary Boron structures? 
First note that our results (Fig. ||) are well fit by a separable form Gz(R, R, 9) = gz(R) 2 Az(9). 
Provided the results are physical reasonable, it is straightforward to interpolate gz and Az 
with respect to R and 9, respectively, and (less easily) with respect to Z. Then one just 
needs to replace the Z definition used in this paper by one of the well-known formulas that 
depends smoothly on the coordinates, e.g. r ij 5 ) 2 / (^2 3 r ij W )- So finally we would set 

Gz{ri,rj,9) = gz(ri)gz(rj)Az(0). in eq. ([!]), well-defined for general structures (provided 
4 < Z < 7 for all atoms); this would reduce to our present results (|^) within each family of 
Z-coordinated uniform structures. 

An attractive application of classical potentials would be in the liquid phase. Here ab 
initio molecular dynamics studies (on a 48 atom system) found detailed results for bond-angle 
distributions and other correlation functions p4| , in good agreement with experiment [ p3[ . 
Those authors also note that Z = 6 but the icosahedra are broken up (the inverted- 
umbrella coordination shell is no longer found), contrary to earlier thought. Potentials would 
offer a chance to model the structure of amorphous a-B, in which few details are known but 
(from radial distribution functions) the B12 icosahedron is believed to be the motif [^ij. It 
would be interesting to explore more ordered (e.g. micro-quasicrystalline) or less ordered 
(like the liquid) variants of this picture. Finally, if our method were extended to extract a 
classical potential from Si and B/Si structures, this could be applied to the icosahedral [^5) or 
cuboctahedral Eq] B12 clusters which are believed to precipitate in B-doped Si. 



( 4 ) The uninverted umbrella was lowest in energy, but there exists no extended structure in which 
all atoms have this coordination shell. 
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To conclude, we have presented a novel approach to generating a database for potential 
fitting which includes so many structures that the angular potential may be fitted rather 
than assumed. We obtain a reasonable description using potential terms that are local to 
the coordination shell of each atom, and which could be extended to general Boron structures. 
Modeling of the speculative or poorly-known quasicrystal, amorphous, unsolved crystal, liquid, 
and Si inclusion forms of Boron would profit from such a potential. 

We thank M. Teter, D. Allen, and J. Charlesworth for providing code and support in the 
LDA calculation, A. Quandt and M. Sadd for comments, K. Shirai and P. Kroll for discussions. 
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of the Cornell Center for Materials Research supported by NSF grant DMR-9632275. 
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